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Rapidity distribution of gluons in the classical field model for heavy ion collisions 
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The rapidity distribution of gluons produced in heavy ion collisions is studied by a numerical 
computation in 2+1-dimensional classical Yang-Mills theory. By assuming that the classical source 
strength g 2 ii depends on rapidity as g 4 fi 2 ~ e ±Xy and studying collisions of two nuclei with different 
g 2 /j, we find that the rapidity distribution of produced gluons at central rapidities is very broad. The 
transverse energy is seen to decrease even more slowly as a function of y than the multiplicity. We 
discuss these results and the range in y and y/s where they are applicable in the light of experimental 
results and other theoretical calculations. 
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INTRODUCTION 



The longitudinal distribution of particles produced in 
a relativistic heavy ion collision is experimentally most 
easily measured as a function of the pseudorapidity 77, be- 
cause measuring pseudorapidity does not require particle 
identification. This measurement has been done for a 
wide range in r\ e.g. by the PHOBOS experiment 0. 
Theoretical calculations, on the other hand, naturally 
produce rapidity distributions. Given the different par- 
ticle species and the different stages of the collision from 
the initial conditions to hadronisation and decoupling, 
finding the right way to transform between pseudorapid- 
ity and rapidity is not always simple. The uncertainty in 
interpreting experimental pseudorapidity distributions in 
terms of theoretical calculations of rapidity dependence 
is discussed in e.g. Q. More recently, however, the 
BRAHMS and STAR Q experiments have measured 
charged particle yields also as a function of the rapidity y. 
Especially the BRAHMS data, covering a wide range in 
y, facilitates the comparison between theoretical models 
and experimental results. 

The aim of this paper is to calculate numerically the 
rapidity dependence of gluon production in the classical 
field model around central rapidities. In Sec.^we briefly 
review this model and how it can be applied to study 
gluon production in heavy ion collisions. In Sec. II I II we 
discuss the relation between the McLerran-Venugopalan 
model and saturation and argue that the rapidity depen- 
dence of particle production can be studied by varying 
the strengths of the classical sources as g 4 [i 2 ~ e ±Ay , at 
least for large enough yfs and for small enough rapidity. 
Numerical results from applying the classical field model 
to collisions of two nuclei with different color charge den- 
sities are presented in Sec. II VI and these results and their 
applicability to heavy ion phenomenology discussed in 
Sec. El 



II. THE CLASSICAL FIELD MODEL 

The McLerran-Venugopalan model for the small x 
wavefunction of an ultrarelativistic nucleus was suggested 
in p| IE ■ The classical field model for the initial stage 
of a collision of two heavy ions, based on the McLerran- 
Venugopalan model for the nuclear wavefunction, was 
formulated in Q and in [ij . 

Let us assume we have two nuclei moving along the 
light cone, corresponding to a current 

J» = 6»+6(x-)p (1) (x T ) + <P-<5(z+)p (2) (x t ). (1) 

The two colour charge densities p( m )(xy) are, indepen- 
dently for the two nuclei, drawn from a random ensem- 
ble, which in the original McLerran-Venugopalan model 
is taken to be Gaussian: 

(pU(*r)p\ m) (yT)) = s 2 M( m) <5 ob £ 2 (xr-yr), m = 1,2, 

(2) 

where fi is a parameter describing the transverse density 
of color charges 1 . 

We then want to find the color fields generated by this 
current using the classical equations of motion 



(3) 



In the light cone gauge (A + = for nucleus (1), A~ = 
for nucleus (2)) one first calculates the pure gauge fields 
corresponding two the two nuclei: 

A\ m) (x T ) = ^ (m) (x T )a i C/ ( t m) ( XT ), m = 1,2. (4) 

These depend on the Wilson lines in the covariant gauge, 

f/(i)(x T ) = Pexpji J dx"yl+ ov (a:",x T )| (5) 

L/ (2 )(xt) = Pexpji J dx + A~ ov (x + ,x. T )^ , 
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1 The relation to the convention introduced in llll is g 2 fi = A s . 
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which can, for infinitely Lorentz-contracted nuclei, be cal- 
culated as 



^(m)( x r) = CX P 



-*5-=S-(*r) 



In a temporal gauge A T = the initial condition at r = 
for the color fields At(t,xt) and A v (t 7 xt) is given by 
these pure gauge fields corresponding to the two nuclei: 



^(0,x T ) = A l cl) (x T )+4 2) (x T ), 
^(0,x T ) = |[4)(x T M{ 2) (x T ) 



the Wilson lines, Eq. J5J, in the adjoint representation 
(see |33| for the detailed calculation). Up to a logarith- 
mic uncertainty which we will neglect in this study, the 
(6) strengths of the classical sources, <? 2 M(i,2)) should thus 
also depend on the x probed in nuclei (1) and (2) respec- 
tively. In the perturbative weak field limit of the model 
we are considering here the partons are produced in a 
2 — » 1 process, with a gluon produced at a given pt and 
y coming from gluons in the initial nuclear wave functions 
at [H: 



(7) 



To find the gauge field in the future light cone r > 
one then solves the gauge field equations of motion using 
these initial conditions. In the gauge A T = it is easy to 
find the Hamiltonian and thus the energy of a given field 
configuration. Additionally, fixing the Coulomb gauge in 
the transverse plane, Vr • At — 0, one can also define 
a gluon multiplicity corresponding to the classical fields. 
The method of solving the classical Yang-Mills equations 
numerically has been developed in 03. First numerical 
results for SU(2) were found in [UllJ and for SU(3) m 
[Til ITsl II?H | 2 The numerical code and notations in this 
paper are those used in |l5j |. 



III. SATURATION AND RAPIDITY 
DEPENDENCE 

The saturation scale Q * is an important concept in 
small x physics (see e.g. [Til IT"9j for a review). It has 
been found to depend on the value of x probed in the 
process as Q 2 ~ x~ x both in analytical studies of the 
quantum evolution of the color sources in the classical 
field model [MISllSlEllllpllpI 3 and in phenomeno- 
logical studies of DIS data |27 l l2a p The value of A found 
in [13 is A = 0.277. . .0.288. Lip*! l3C| the values used 
are A = 0.25 . . .0.3. In presenting the numerical results 
in Sec. IIVI we shall keep A arbitrary and use it only in 
Sec. ]V\ when comparing with experimental results and 
other theoretical calculations. 

In the McLerran-Venugopalan model the saturation 
scale Q s is related to the strength of the classical sources 
by 



Ql 



Ait 



In 



ffV 

A QCD, 



(8) 



Note that this is the saturation scale as defined in e.g. 
[Tfj| , and differs from the one used in the context of DIS 
(e.g. 26]) by the color factor Ca vs. Cf 4 - The rela- 
tion JSJl can be derived by calculating the correlator of 



2 See also the erratum Il7l . 

3 See e.g. for a more comprehensive list of refe rences. 

4 This difference is discussed in more detail in e.g. Illll.'jll . 



Xi2 = Pre 



±v 



(9) 



Assuming that the dominant transverse momentum scale 
(pt) does not depend on y one is lead to assume that to 
calculate gluon production at a rapidity y one must take 



/*(!) 



= fi 2 e Xy and 



M(2) = M 



-Ay 



(10) 



where M(i,2) are the source strengths appearing in the 
correlator in Eq. J2J. Note that the geometric mean 
/i = y/n(i)p.{2) is the quantity appearing in Eqs. illlfl . 
(fT2"|) . Our assumption that (pr) does not depend on y is 
valid only at small y in AA-collisions where, by symme- 



try, one expects (pt) 



;i + 0(AV)) (PT) y =o- Further 



away from central rapidities one could assume (pr) to 
depend on one of the saturation scales of the two nuclei, 
and the relation (|10f) cannot be used any more. A sepa- 
rate question from the selfconsistency of the calculation 
discussed above is whether y/s at RHIC is high enough 
for this model to be applicable; we will return to this 
question in Sec. 



IV. RESULTS 

We calculate the transverse energy and multiplicity per 
unit rapidity as in |15|. and show our results in terms of 
the dimensionless ratios 



ft 



1 



AE 



g A H 3 TrR 2 A dri 



and 



dN 



N 



g 2 IJ?%R 2 A dr] 



(11) 



(12) 



These depend on Xy through the source strength (by 
Eq. iJTUJl) and the dimensionless parameter characteris- 
ing the field strength g 4 fj, 2 irR A . We shall consider g = 2 
and ttR a = 140 fm 2 as constants and vary the source 
strength by using different values of /i. As in [l^, the 
numerical computation presented here is done assuming 
cubic nuclei, whose projection to the transverse plane fills 
the whole 2-dimensional lattice (of transverse area ttR a ), 
enabling the use of periodic boundary conditions. 

Figure n shows the behaviour of Jn and Fig. [21 of Je as 
a function of Ay for transverse lattices of 256 2 and 512 2 
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Multiplicity (f ) for different field strengths 



Multiplicity (f ) for different field strengths 



256 -lattice 



o |i = 0.3GeV 
■ |i = 0.5GeV 
♦ |i = 0.8 GeV 
-0.3 GeV, fit 
0.5 GeV, fit 
- 0.8 GeV, fit 




512 -lattice 

~i 1 I 



o |i = 0.3GeV 
■ |l = 0.5GeV 
♦ |i = 0.8 GeV 
-0.3 GeV, fit 
0.5 GeV, fit 
- 0.8 GeV, fit 




Figure 1: The dimensionless ratio /jv for different rapidities on two different lattice sizes (left: 256 2 , right: 512 2 ) and for 
different field strengths. Also shown are the Gaussian fits with the widths in Table [I] 
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Figure 2: The dimensionless ratio /s on two different lattice sizes (left: 256 2 , right: 5f2 2 ) and for different field strengths. 



points. One can see that whereas the multiplicity de- 
creases slowly with rapidity, the transverse energy does 
not change significantly, meaning that the energy per par- 
ticle increases. This increase is so slow, however, that our 
approximation in Sec. 11111 of considering (px) roughly in- 
dependent of y is justified for Xy < 1. The behavior of 
the gluon spectrum is illustrated in Fig. |3 which shows 
the ky-distibution of gluons for two different rapidities. 
This hardening of the spectum is an inherent feature in 
a distribution depending on two transverse momentum 
scales gV(i,2)- Note that the presence of two widely dif- 
ferent momentum scales also poses difficulties for the nu- 
merical calculation, as the harder scale approaches the 
lattice ultraviolet cutoff and the softer scale approaches 
the infrared divergent weak field result. For this reason, 
although one can use this method to study gluon spec- 
tra and the Cronin effect in "pA" -collisions (with one 



of the sources very weak), as has been done in |l6j |. it 
is perhaps not realistic to calculate integrated multiplic- 
ities for the "pA"-case. Instead, saturation physics in 
"pA" -collisions 5 can be studied using a qualitatively dif- 
ferent treatment for the proton and the nucleus (see e.g. 

00 El S3) 

We fit the rapidity dependence of the multiplicity with 
a Gaussian form 

fN = /ocxp (-J^) = /o exp (-|^) • (13) 

Because the result only depends on the combination Ay, 



5 Actually the relevant RHIC experiments, see e.g. 
are d+Au-collisions. 
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1 2 3 4 5 6 
kJg ii 

Figure 3: The spectra dN/ d 2 kr of produced gluons for 
a symmetric collision = M(2)) an d an asymmetric one 

(fJ>n) < M(2)) plotted as a function of k/g 2 /i, where k 2 — 
4^ sin 2 kill. 







Act, 256 2 -lattice 


Act, 512 2 -lattice 


0.3 GeV 


5184 


1.99 ± 0.02 


1.96 ± 0.03 


0.5 GeV 


14400 


1.56 ± 0.02 


1.66 ± 0.02 


0.8 GeV 


36864 


1.32 ± 0.03 


1.35 ± 0.02 



Table I: The widths of the Gaussians Act fitted to the multi- 
plicities in Fig. 



with A arbitrary so far, the result from the fit is actually 
the combination Act. The values of Act from the fits are 
listed in Table [I] One sees that there is a small depen- 
dence lattice size and a surprisingly strong dependence 
on the field strength g 4 /j, 2 TrR\. The width Act seems to 
depend on the field strength more than = 0) and 

/jv(y = 0). This means that Act is more sensitive to the 
lattice infrared cutoff than the energy and multiplicity 
themselves. Note that at least in this context there is no 
fundamental reason for the Gaussian form, Eq. (| 1 3|) : by 
symmetry the distribution must be even in y and here we 
are mainly interested in the second derivative at y = 0. 
Because of the numerical reasons discussed above we do 
not want to go to very large values of y and a form like 
like [1 + y 2 / (2na 2 )]~ n for some n would be just as good. 
We choose the Gaussian because it has also been used in 
experimental Q and other theoretical studies [4^. The 
rapidity distribution of the energy does not seem to nat- 
urally lend itself to a simple fit that would illustrate its 
structure any more than Fig. [21 



V. COMPARISON WITH EXPERIMENTAL 
RESULTS AND OTHER THEORETICAL 
CALCULATIONS 

It was found in 15] that assuming ideal hydrody- 
namical expansion and consequently entropy conserva- 
tion RHIC multiplicities are best reproduced by taking 
g 2 /i 2 GeV i.e. fj, « 0.5 GeV. Taking the correspond- 
ing value Act ss 1.66 from TableQJand A ~ 0.25 . . . 0.3 one 
gets cr rs 5.5 . . . 6.6. 

The dependence of the saturation scale on rapidity has 
been exploited to study heavy ion phenomenology also by 
Kharzeev & Levin 29]. The authors use a saturation- 
inspired gluon distribution and perform a pertubativc 
calculation of gluon production. The same approach has 
also been used by Hirano & Nara [3(j to obtain initial 
conditions for a hydrodynamical calculation. Estimating 
from Fig. la of [3(| the width of the distribution seems 
to be ct ss 3, but the shape seems to be flatter than a 
Gaussian of this width at small y and fall more rapidly 
at large y than a Gaussian. 

Our calculation differs from these papers in that we are 
performing a numerical calculation by soving the classi- 
cal field equations, not a perturbative calculation using 
unintegrated gluon distributions. One advantage in our 
calculation is that the multiplicity comes out naturally 
as a smooth function in y around y — 0, without the 
e~ A ! y !-discontinuity in the first derivative in 29]. In a 
perturbative calculation one can try to incorporate fea- 
tures of high-cc physics, such as the (1 — x) 4 -behaviour 
of the parton distributions imposed by hand in |29l l30j , 
whereas our calculation stays within the framework of a 
small x model without including this kind of effects. This 
means that our calculation is limited to regions around 
y = 0, where the fields of both nuclei are strong and can 
be treated classically, and such values of that large x 
effects are not important at midrapidity. 

These results can also be compared to the 
pQCD+saturation model calculation 0] result of 
ct = 5.9. Whereas in 0] a increases slightly when the 
saturation scale grows, in our calculation it decreases 
and the distribution becomes more peaked. 

It is not straightforward to compare this calculated 
initial state gluon distribution to measured rapidity dis- 
tributions in the final state. The approach of |29j is 
to assume that they are equal and only correct for the 
transformation between rapidity and pseudorapidity; in 
the hydrodynamical calculation of [43| the distribution 
broadens, but only slightly, during the hydrodynamical 
evolution. The result of the BRAHMS collaboration 
for the rapidity distribution of charged pions is ct « 2.3, 
which is close to the result of [3(j and considerably less 
than our result. This could be interpreted as an indica- 
tion that the main features of the rapidity dependence 
of particle production in AA-collisions at RHIC are not 
dominated by saturation physics but by the large x be- 
havior of the parton distribution functions that cause 
the multiplicity to fall faster for large rapidities than the 
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Figure 4: The first dashed curve labeled "u = 6" is the Gaussian e y // ' 2 ' 6 The second curve, labeled "(1 — x) 4 , a = 6" is 
the same Gaussian multiplied by (1 - e M /200) 4 (l - e" y /200) 4 , where the factor 1/200 comes from assuming y/s/A = 200 GeV 
and (pt) ~ 1 GeV. The solid curve is e~ v ^ 2 ' 3 a Gaussian with width a = 3. This figure is not intended as a reproduction 
of the calculation of |2!JI. but rather as a qualitative demonstration. 



classical field model would suggest. The crucial impor- 
tance of the factor (1 — x) 4 in the distribution functions 
of |29, y(| is illustrated in Fig. 0] which shows how, at 
RHIC energies, this factor can transform a broad Gaus- 
sian with a = 6 into a form that is closer to a = 3. 
Note that the situation at the LHC will be quite differ- 
ent as the values of x probed at central rapidities will be 
an order of magnitude smaller than at RHIC and thus 
much less influenced by large x physics. Another way of 
viewing this is that energy conservation limits the size 
of a rapidity plateau at central rapidities and forces the 
multiplicity to decrease faster with increasing \y\ 

In this calculation the average transverse momentum 
of the gluons, (px) — ^g 2 ^ increases for larger rapidi- 
ties. The experimentally measured mean px in the final 
state is approximately constant or decreases. Although 
one can, in the hydrodynamic scenario, argue that a large 
part of the transverse energy goes into the longitudinal 
expansion of the system, it is hard to understand how 
the energy could decrease by a larger amount for larger 
rapidities. In the calculation of |3(j the transverse energy 
decreases faster for larger y due to the explicit (1 — x) A - 
factor and momentum cutoffs in the transverse integra- 
tion. 



VI. CONCLUSION AND OUTLOOK 

We have presented a numerical calculation of the ra- 
pidity dependence of gluon production in heavy ion colli- 
sions. The calculation is performed in a 2 + 1-dimensional 
classical Yang-Mills model assuming that the strength of 
the classical sources vary with rapidity according to the 
simple relation g 2 fj, ~ e ±Xy . Our result is that the multi- 



plicity is approximately a very broad Gaussian in rapid- 
ity with a width a « 6. The distribution resulting from 
the classical field model is broader than the one observed 
experimentally. This could indicate that the rapidity dis- 
tribution of particle production in heavy ion collisions at 
RHIC experiments depends more on the (1— x) -like large 
x behavior of parton distributions than actual saturation 
physics, i.e. the Q\ ~ ir~ A -dependence of the saturation 
scale. At the LHC the situation could be quite differ- 
ent. Note, however, that our calculation only applies to 
central rapidities in AA-collisions and does not address 
pA-collisions. 

The transverse energy produced in this model does not 
decrease for larger rapidities as fast as one would physi- 
cally expect. This can be due to the model not including 
some relevant large x physics. Another potential reason 
could be that a 2+1-dimensional model is not sufficient to 
capture all the aspects of the longitudinal dependence of 
the problem, as is the case e.g. for quark pair production 
p4| . This will be understood better if one is able to per- 
form a full 3+1-dimensional classical field computation 
with a more detailed inclusion of the quantum evolution 
of the sources. This calculation would also be important 
for understanding the subsequent thermalisation of the 
gluons produced in the collision. 
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